package edu.gmu.atelier;

/**
 * Performs SVD factorization on the passed matrix. 
 * SVD:  AV   =   U?  where V and U are orthonormal
 *       A[v0,v1,..vr...vn] = [u0,u1,...ur...un][?1,?2,...?r,0,0,...0n]
 *       A = U?V^T
 * 
 * To determine orthonormal matrix V (make U's go away)
 * A^T A =(V ?^T U^T)(U ? V^T)= V ?^2 V^T
 * 
 * To determine orthonormal matrix U (make V's go away)
 * A A^T =(U ? V^T)(V ?^T U^T)= U ?^2 U^T
 * 
 * @author James H. Pope
 */
public class SVDFactor
{
    private Matrix a  = null;
    private Matrix u  = null;
    private Matrix s  = null;
    private Matrix vt = null;
    
    public SVDFactor( Matrix a )
    {
        this(a, QRFactorFactory.HOUSEHOLDER, EigenQRFactor.DEFAULT_MAX_ITERATIONS);
    }
    
    public SVDFactor( Matrix a, QRFactorFactory f, int maxIterations )
    {
        this.a = a;
        
        Matrix ata = a.transpose().mult(a);
        
        EigenQRFactor fv = new EigenQRFactor(ata, f, maxIterations);
        //System.out.println("Iterations = "+fv.getIterations() + " converged "+fv.converged());
        //System.out.println("A^T A: ?^2 = "+fv.getD());
        //System.out.println("A^T A: V   = "+fv.getS());
        Matrix v = fv.getS();
        
        Matrix aat = a.mult(a.transpose());
        EigenQRFactor fu = new EigenQRFactor(aat, f, maxIterations);
        //System.out.println("Iterations = "+fu.getIterations() + " converged "+fu.converged());
        //System.out.println("A A^T: ?^2 = "+fu.getD());
        //System.out.println("A A^T: U   = "+fu.getS());
        this.u = fu.getS();
        
        // Sanity check is to make sure both derived ?^2 are the same
        Matrix d = fu.getD(); // may not be "perfectly" diagonal
        // Go through make "close" to zero entries exactly zero so can take sqrt
        Epsilon e = Epsilon.makeDecimal(14);
        for (int i = 0; i < d.getVectors().size(); i++)
        {
            Vector vector = d.getVector(i);
            for (int j = 0; j < vector.getLength(); j++)
            {
                Element elem = vector.getElement(j);
                Element zero = elem.getZero();
                if( elem.similar(zero, e) )
                {
                    vector.setElement(j, zero);
                }
            }
        }
        //System.out.println("D="+fu.getD());
        /*
         * If not perfectly diagonal, this will exception, fine, what about
         * negative numbers!  Would seem so work with complex numbers.
         * Need to confirm - negative eigen values in sigma
         * require complex numbers or is valid exception? Said another
         * way, can A A^T produce a matrix with non-real numbers?
         */
        this.s = d.sqrt(); 
        
        // Not sure if we should keep v or vt or both?
        this.vt = v.transpose();
    }
    
    public Matrix getA()
    {
        return this.a;
    }
    
    public Matrix getU()
    {
        return u;
    }
    
    public Matrix getS()
    {
        return s;
    }
    
    public Matrix getVT()
    {
        return vt;
    }
}
